Finite-temperature Simulations for Magnetic Nanostructures 
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We examine different models and methods for studying finite-tempera-ture magnetic hysteresis 
in nanoparticles and ultrathin films. This includes micromagnetic results for the hysteresis of a 
single magnetic nanoparticle which is misaligned with respect to the magnetic field. We present 
results from both a representation of the particle as a one-dimensional array of magnetic rotors, and 
from full micromagnetic simulations. The results are compared with the Stoner-Wohlfarth model. 
Results of kinetic Monte Carlo simulations of ultrathin films are also presented. In addition, we 
discuss other topics of current interest in the modeling of magnetic hysteresis in nanostructures, 
including kinetic Monte Carlo simulations of dynamic phase transitions and First-Order Reversal 
Curves. 
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I. INTRODUCTION 



Although hysteresis in fine magnetic particles has been 
intensively studied for many years, there is currently sig- 
nificant interest in reexamining our understanding of this 
phenomenon. Partly, this interest is driven by the po- 
tential application of hysteresis in nanostructures to new 
technologies such as Magnetic Random Access Memory 
(MRAM) and ultra-high-density magnetic recording. For 
the past several years, the areal density of hard drives has 
been doubling every 18 months, and is rapidly approach- 
ing the limits of conventional longitudinal recording tech- 
nology. At the same time, data rates in these drives have 
increased significantly, with the serial interface standard 
at the time of writing providing a peak data transfer 
rate of 2.4 Gb/s^. This has led the magnetic record- 
ing industry to look at new recording paradigms such as 
patterned media and self-assembled arrays of nanostruc- 
tures. In fact, the first laptop computer incorporating a 
hard drive based on perpendicular recording technology 
was recently introduced^. It is crucial, then, to under- 
stand the complex process of hysteresis in these systems. 

At the same time, recent advances in computational 
ability, both in terms of new algorithms and available 
computer resources, allow for numerical studies never be- 
fore possible. Plumer and van Ek 3 , for instance, have 
studied the effects of anisotropy distributions in perpen- 
dicular media using a micromagnetic model. Their re- 
sults (Fig. 1) show how anisotropy distributions tend 
to reduce the squareness of the loop and, therefore, the 
signal to noise ratio (SNR). Gao et al. have recently car- 
ried out similar studies of tilted perpendicular media 4 
and polycrystalline media™. Another important effect 



which can be better understood through simulations is 
Barkhausen noise. This effect also decreases SNR, partic- 
ularly in new thin-film media with soft underlayers. Dah- 
men, Sethna, and coworkers used a random-field Ising 
model to examine the origins of Barkhausen noise and 
have been able to relate it to avalanches and disorder- 
induced critical behaviosSi^. These results illustrate two 
of the ways simulations can be used to help understand 
both fundamental processes in hysteresis and their appli- 
cations to new technology. 

Here we present an overview of several common ap- 
proaches to studying hysteresis in magnetic nanostruc- 
tures. We then present results of large-scale computer 
simulations of hysteresis in single iron nanoparticles when 
the magnetic field is misoriented with respect to the long 
(easy) axis of the elongated particles. We also examine 
other recent advances in the study of magnetic hystere- 
sis, such as kinetic Monte Carlo simulations of dynamic 
phase transitions and First-Order Reversal Curves. 



II. MODELS AND METHODS 

A. Coherent Rotation 

Given a single-domain particle with uniaxial 
anisotropy, it is possible to find the metastable 
and stable energy positions of the magnetization when a 
magnetic field is applied at an angle to the easy axis. It 
is assumed that the magnetization can be represented 
by a single vector M, with constant amplitude, Ms- 
The energy density of the system is then 



E = K sin 2 d - M S H cos(</> - 0) , 



(1) 
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where K is the uniaxial anisotropy constant, H is the 
magnetic field applied at an angle <fi to the easy axis, and 
•d is the angle the magnetization makes with the easy 
axis. Stoner and Wohlfarth showed that for coherent 
reversal of the magnetization, the spinodal curve beyond 
which the metastable energy minimum disappears and 
switching occurs is given bjal, 
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where Hax and Hay are the respective components of the 
magnetic field normalized by the anisotropy field Hk = 
2K/M$, along the easy and hard axes. Equation (2) is 
the well-known equation of a hypocycloid of four cusps, 
also known as an astroid. 



B. Micromagnetics 

For systems in which the spins are not aligned and/or 
the field is changing too rapidly for the magnetization 
to reach its quasi-static value, it is usually necessary to 
use a non-perturbative technique such as micromagnet- 
ics to describe the reversal process. The basic approach 
is to divide the system into a coarse-grained set of sites. 
Each site is associated with a position f , and its magne- 
tization is represented by a single magnetization vector 
M (f), whose norm is the saturation magnetization Mg, 
corresponding to the bulk material (a valid assumption 
for temperatures well below the Curie temperatureiS). 
The time evolution of each spin is given by the Landau- 
Lifshitz-Gilbert (LLG) equation 11 ! 12 . 13 , 



dM(r-) 
dt 



7o M(fi) x 



H(n) 



7o Ms 



dM(r-) 
dt 



(3) 

where H (r*,) is the total local field at the i-th site, 70 is 
the gyromagnetic ratio (1.76 x 10 7 rad/Oes), and a is a 
dimensionless damping parameter which determines the 
rate of energy dissipation in the system. The first term 
represents the precession of each spin around the local 
field, while the second term drives the magnetization to 
align with the field. The LLG equation can easily be 
rewritten in a form more convenient for numerical into- 
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(4) 



For the sign of the undamped precession term, we follow 
the convention of Brownii. 

The total local field, H (fi), controls the dynamics and 
contains all of the interactions between each site and the 
rest of the system; it is defined by 

H(n) = — P^. (5) 
dM(f) 



Here, Ei is the free energy of the i-th site and the op- 
erator d/dM{fi) = (d/dM x (fi))x+ (d/dM y {f t ))y + 
(d/dM z (fi)) z. The different terms that contribute to 
H (fi) combine via linear superposition, 

H(fi) = H z (f)+H e (f)+H D (f l )+H a (f l )+H n (f) . (6) 

Here, Hz (f) is the externally applied field (Zeeman 
term), H e (f) is due to exchange interactions, Hp (f) 
is the dipole field, H a (f) is the anisotropy field (in our 
simulations taken to be zero), and H n (f) is a random 
field representing the effects of thermal noise. 

The exchange contribution to the local field repre- 
sents local variations between the alignment of M (f) 
and neighboring sites and can be represented by 
l\ V 2 M (r*i)i^. In our simulations, this is implemented 
by 
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where the summation is over the nearest neighbors of 
f, rii is the number of neighbors of site i, and the term 
riiM (fi) is included so that H e = when all of the spins 
are aligned. The exchange length, l e , is defined in terms 
of the exchange energjii^, E e — - (^/2) / drM ■ V 2 M, 
in a continuous system. For our discrete system of finite- 
sized cells, this means the magnetization can be viewed 
as rotating continuously from the center of one cell to the 
center of each neighboring cell along the line joining the 
two. 

At non-zero temperatures, thermal fluctuations con- 
tribute a term to the local field in the form of a stochas- 
tic field H n (f), which is assumed to fluctuate indepen- 
dently for each spin. The fluctuations are assumed to be 
Gaussian, with zero mean and (co)variance given by the 
fluctuation-dissipation theorer^- 1 1 ' 



(H^(f,t)H n ^(f t ,f)) = ^^S(t-t')6^S iti > , (8) 

7 iWs V 

where H n ^ (f) indicates one of the Cartesian coordinates 
of H n (f). Here, T is the absolute temperature, fcs is 
Boltzmann's constant, V = (Ar) 3 is the discretization 
volume of the numerical integration, and 8 r y is the Kro- 
necker delta representing the orthogonality of the Carte- 
sian components. Although this result was derived for 
an isolated particle, recent work by Chubykalo, et al. 
indicates that this result will hold for interacting sys- 
tems as welli£. In this paper, we present micromagnetic 
results for two different models. The first model is of 
a nanoparticle with dimensions 5.2nmx5.2nmx88.4nm. 
The cross-sectional dimensions are small enough (« 2 
l e ) that the assumption is made that the only sig- 
nificant inhomogeneities occur along the long axis (z- 
direction)i2ii£. The particles of this model are therefore 
discretized into a linear chain of 17 spins along the long 




FIG. 1: The effect of anisotropy distributions on hysteresis in perpendicular media. A Gaussian distribution of anisotropy fields 
is used with a mean value of H^u = 7000 Oe. In the bottom left corner, the curves correspond to standard deviations in Hk 
of (from left to right) 2000 Oe, 1000 Oe, 500 Oe, and OOe. The saturation magnetization is 350 emu/cm 3 . Data courtesy of M. 
Plumeu- 



axis of the particle. We will call this model the stack-of- 
spins model. 

In this simple model, the local field due to dipole-dipolc 
interactions is calculated asi£*i£ 

g (^ = (Ar)3^ 3 ^-^))-^) , (9) 

where is the displacement vector from the center of 
cube i to the center of cube j, and is the correspond- 
ing unit vector. The volume factor (Ar) 3 results from 
integration over the constant magnetization density in 
each cell. The second model, which we will refer to as the 
full micromagnetic model, simulates a single nanoparticle 
with dimensions 9 nm x 9 nm x 150 nm. The dimensions 
were chosen to correspond to arrays of nanoparticles fab- 
ricated by Wirth, et ali&. In this model, the system is 
discretized into 4949 sites (7 x 7 x 101) on the computa- 
tional lattice. The size of the system makes calculation of 
dipole interactions in the conventional manner (as done 
for the stack-of-spins model) computationally impracti- 
cal. It is therefore necessary to use a more advanced 
algorithm to make the simulation tractable. 

The two most popular choices are the traditional Fast 
Fourier Transform (FFT) and the Fast Multipole Method 
(FMM)i2u Here, we used the Fast Multipole Method, the 
exact implementation of which is discussed elsewhere^, 
because it has several advantages over the FFT. The 



biggest difference is that the FMM makes no assump- 
tions about the shape of the underlying lattice, while 
the FFT assumes a cubic lattice with periodic bound- 
ary conditions. The consequence of this is that numer- 
ical models of systems without periodic boundary con- 
ditions which use the FFT require empty space around 
the system so that the boundary conditions do not affect 
the calculation. The FMM requires no such "padding" . 
Furthermore, the FFT requires O (NlnN) operations to 
calculate the magnetic scalar potential (from which the 
dipole field is calculated). The FMM algorithm, while it 
has a larger computational overhead, requires only O (N) 
operations for the same calculation. This means that, 
while the FFT is a good choice for small cubic lattices, 
the FMM is better for large, incomplete, or irregular lat- 
tices. The public-domain psi-Mag toolset now provides a 
flexible implementation of the FMM designed for use on 
high performance, parallel computersSS. 

Material properties in both models were chosen to cor- 
respond to bulk Fe. The saturation magnetization is 
1700 emu/cm 3 (kA/m) and l e = 2.6 nm. We take the 
damping parameter a = 0.1 to correspond to the un- 
derdamped behavior usually assumed to be present in 
nanoscalc magnets. Although this value is approximately 
an order of magnitude larger than the value obtained ex- 
perimentally using Ferromagnetic Resonance (FMR), it 
has been noted that the FMR value is for small devia- 
tions of the magnetization from equilibrium and is not 
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representative of the large deviations which occur dur- 
ing reversal— In general, care should be taken in es- 
tablishing an appropriate damping parameter to use in 
simulating a particular magnetic nanostructure, as also 
illustrated in other recent micromagnetic studies^ 2 * 2 ^*— 
It is worth noting that, even in systems which reverse 
coherently at high speed, deviations from the quasi-static 
Stoner-Wohlfarth (SW) model will be expected. He, et 
al— i showed that for square-pulse fields with fast rise 
times (< 10 ns) and small values of the damping con- 
stant (< 0.2), the shape of the astroid changes. The 
result is that the minimum switching field is reduced be- 
low the SW limit of 0.5 Hk , and the angular dependence 
is no longer symmetric around 45°. For the frequencies 
and damping parameter used here, the deviations from 
the SW model are small (< 5 percent difference in the 
switching fields) and may be neglected. 



C. Monte Carlo Simulations: Kinetic Ising and 
Heisenberg Models 

A second approach to modeling the dynamics of mag- 
netic systems involves Monte Carlo techniques, which 
have been applied to a wide variety of systems since 
their introduction by Metropolis, et al. 26 . As described 
above, the micromagnetics approach uses a stochastic 
(i.e. random number-based) variable to introduce ran- 
dom fluctuations into an otherwise deterministic system. 
In contrast, Monte Carlo simulations are fully stochastic 
and proceed by considering possible transitions between 
states of the system and executing these transitions with 
a probability which depends on the system's energy and 
temperature. 

The static Monte Carlo algorithm consists of a re- 
peated three-step process. First, choose a (pseudo- 
random number. The random numbers chosen may be 
uniformly distributed, or they may be chosen based on 
a particular probability distribution that depends on the 
specific simulation to be performed. Second, choose a 
trial move from the current state to a new state. Third, 
accept or reject the trial move depending on the ran- 
dom number and some acceptance rule consistent with 
the problem under consideration. 

Consider the ferromagnetic Ising model on a regular 
lattice with periodic boundary conditions. Each site on 
the lattice has a spin which can align either parallel or 
anti-parallel to the applied field and takes on values of 
Si = ±1 accordingly. The energy of the Ising lattice is 
then 



E = -jJ2SiSj-H(t)J2Si 



(10) 



where the exchange constant J > is in units of energy, 
and H (£) is the externally applied, time-dependent mag- 
netic field (which in Monte Carlo simulations is custom- 
arily given in units of energy, thus absorbing the mag- 
netic moment per site, jj). The first sum in (10) is over 



nearest-neighbor pairs, while the second sum is over all 
spins on the lattice. 

The static Monte Carlo procedure described above al- 
lows the calculation of equilibrium quantities such as the 
internal energy, susceptibility, specific heat, and magne- 
tization. In order for the lattice to explore each of its pos- 
sible states with probabilities corresponding to the equi- 
librium thermal distribution, the acceptance rule chosen 
must satisfy the condition of detailed balance 27 . Two 
common choices are the Metropolis^ and Glauber— ac- 
ceptance rules. Note that near the critical point, com- 
putation with these simple acceptance rules slows down 
dramatically, and it is therefore useful to use more ad- 
vanced algorithms (such as cluster algorithms 27 ) to cal- 
culate equilibrium quantities. 

In equilibrium calculations, no physical interpretation 
is ascribed to the intermediate spin flips. If, instead, we 
consider the individual spin flips as representing physi- 
cal fluctuations due to the interactions between the spins 
and a heat bath, then the underlying transitions model 
the actual dynamics of the system and acquire a physical 
significance. This application of Monte Carlo simulations 
is known as kinetic Monte Carlo. The random nature of 
the events due to the interaction of spins dictates that 
the spin to attempt to flip must be chosen at random. 
In this paper, we use the Glauber^ acceptance rule, ac- 
cording to which each attempted spin flip is accepted 
with probability 



W = 



exp(-/3Ag t ) 
1 + exp (-/3A^ 



(11) 



Here, AEi is the change in energy that results if the pro- 
posed flip of the i-th spin is accepted, and (3 — (ksT) 
With a uniformly distributed random number, r G [0, 1], 
a randomly chosen spin is flipped if r < W. Each po- 
tential spin flip is considered a Monte Carlo step. The 
basic time step of the Monte Carlo process is measured in 
Monte Carlo Steps per Site (MCSS). This time is related 
to the algorithm and in general is only approximately 
proportional to the physical time of the system. Recently, 
however, progress has been made in connecting analyti- 
cally the MC simulation time to the simulation time of 
the Langevin-based micromagnetic techniques discussed 
above, for which there is a clear relationship to physical 

time 29,30,31,32 _ 

The Glauber dynamic of (11) can be derived from a 
quantum spin-i Hamiltonian coupled to a fermionic heat 
bath— Recently, other dynamics have been derived from 
coupling a quantum spin-i system to a phonon heat 
bath— Note that in kinetic Monte Carlo calculations, 
algorithms (such as the cluster algorithm) that change 
the underlying dynamic cannot be used. However, ad- 
vanced algorithms that achieve very large speedups while 
remaining true to the underlying dynamics are possi- 
ble^. It has recently been shown that physically relevant 
functional forms for W can lead to dramatically differ- 
ent values of dynamical quantities such as lifetimes of 
metastable states— 
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It is important to realize that the Monte Carlo tech- 
niques can be applied to other systems as well. Unlike 
the Ising model, the Heisenberg model allows the spins 
to assume any angle with respect to neighboring spins 
and the applied field. The energy of a regular lattice of 
Heisenberg spins with periodic boundary conditions is 

E = — J y^(Si x Sj X + SiySjy + Si z Sj z )—H y ^Sj cos(flj) , 

(id) > 

(12) 

where Si X , Si y , and S% z are the Cartesian coordinates of 
the vector spin Si (with magnitude Si = 1), and 9i is 
the angle between the applied field, H, and the z-th spin. 
As in the Ising model, the first sum is taken over nearest 
neighbors and represents the exchange interactions, while 
the second sum is taken over all spins in the system, and 
represents the interactions of the spins with an externally 
applied magnetic field (Zeeman energy). 

The dynamic consists of randomly choosing a spin to 
update, randomly choosing a new spin direction (either 
uniformly distributed over the sphere or over a cone near 
the current spin direction), and using a Metropolis or a 
heat-bath rate to decide whether to effect a transition 
to the new spin direction. The rate depends on the en- 
ergy difference between the spin configurations as in, for 
example, (11). 

In kinetic Monte Carlo, it is possible to implement the 
algorithm in a rejection-free manner, so that every al- 
gorithmic step performs an update. In this case, each 
algorithmic step in general advances the system by a dif- 
ferent amount of time. For models, such as the Ising 
model, with discrete state spaces this is called the n-fold 
way algorithm— It is possible to make a precise con- 
nection between these the n-fold way and the standard 
implementation of kinetic Monte Carlo— Recently such 
rejection-free methods have been implemented for mod- 
els with continuous state spaces, such as the Heisenberg 
model—, and the efficiency of rejection-free methods in 
various systems has been studied— 



III. RESULTS OF MICROMAGNETIC 
SIMULATIONS 

In this section we summarize recent simulation results 
for magnetization reversal in iron nanopillarssi, and fur- 
ther evaluate these results in light of additional experi- 
mental data on such reversal. 

Figure [2t shows hysteresis loops at T — 100 K for the 
full micromagnetic model with the field misaligned at 0°, 
45°, and 90° to the long axis of the particle. The loops 
were calculated using a sinusoidal field with a period of 
25 ns, which started at a maximum value of 10,000 Oe 
(800 kA/m). In all the loops in this section, the reported 
magnetization is the component along the long axis (z- 
axis) of the particle. Simulations for the full micromag- 
netic model were performed over one half of the period 
and the results reflected to give the full hysteresis loop. 



Consider the case with the field and particle aligned 
(0°). Initially, the large magnetic field tends to align the 
spins with the easy axis. As the field is decreased, the 
spins relax, and the magnetization decreases by approx- 
imately 2%. Eventually, reversal initiates at the ends as 
previously reported— 

At 45° misalignment between the particle and the field, 
the magnetization is initially pulled away from the long 
(easy) axis by the large magnetic field. As the field is 
swept toward zero, the magnetization relaxes until it 
essentially reaches a maximum value of approximately 
0.91Ms at zero applied field. Thermal fluctuations along 
the length of the particle prevent the magnetization from 
reaching saturation. As in the case of 0° , reversal again 
begins by nucleation at the ends of the particle, with the 
growth of these nucleated regions leading to the reversal 
of the particle. Figure |3t shows the z-component of the 
magnetization at selected times during the reversal pro- 
cess for the 45° hysteresis loop of Fig.|2^. It is important 
to note that the particles do not have a uniform magneti- 
zation during the reversal process, even though they are 
single-domain particles. 

For 90° misalignment, the reversal mechanism is quite 
different. The hysteresis loop in Fig. |2K shows that the 
magnetization is essentially perpendicular to the easy di- 
rection until the field reaches a particular value. As the 
field is decreased further, the magnetization relaxes to- 
ward the easy axis. Since nothing breaks the up/down 
symmetry of the system when the applied field has no 
component along the easy axis, the relaxed magnetiza- 
tion can be directed toward either the positive or negative 
z-axis. Figure|3p shows the z-component of the magneti- 
zation for the 90° misalignment at selected times during 
the hysteresis loop of the full micromagnetic model. For 
this case, the nucleation occurs along the entire length of 
the particle, except at the ends. The large demagnetiz- 
ing fields present at the ends (involved in nucleation at 
smaller angles) retard relaxation along the easy axis. 

The hysteresis loops for the stack-of-spins model, 
shown in Fig. [5p, are qualitatively similar to those of 
Fig. Loops at 0°, 75°, and 90° misalignment are 
shown. There are important differences between the two 
models, however. First, without lateral resolution of the 
magnetization across the cross-section, these particles 
exhibit ringing due to the precessional dynamics. Evi- 
dently, the precession of individual moments in the full 
micromagnetic model does not lead to precession of the 
end-cap moment; possibly the spin waves rapidly damp 
out the gyromagnetic motion. 

A second, and more prominent, difference between the 
models is observed in the angular dependence of the 
switching field, H sw , shown in Fig. 0] Here, H sw is de- 
fined as the applied field at which M z is reduced to 0. 
The stack-of-spins model (circles) shows a shape qualita- 
tively similar to what is expected from Stoner-Wohlfarth 
(SW) theory, with a minimum H sw near 45°. The dashed 
curve is the SW theory with Hk = 1600 Oe (for com- 
parison reasons Hk was chosen to be much smaller than 
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(a) Hysteresis loops for the full micromagnetic 
model at T = 100 K with a sinusoidal field of period 
25 ns and maximum applied field of 10 kOe 



(b) Hysteresis loops for the stack-of-spins model at 
T = 10 K with a sinusoidal field of period 200 ns and 
a maximum applied field of 5 kOe 



FIG. 2: Hysteresis loops for (a) the full micromagnetic model with 0°, 45°, and 90° misalignment and (b) the stack-of-spins 
model with 0°, 75°, and 90° between the applied field and the long axis of the particle 



the 10 4 Oe expected for these particles assuming SW be- 
havior). The full micromagnetic model (diamonds), on 
the other hand, has its minimum H sw at 0°, and H sw 
increases as the misalignment angle is increased. Fig- 
ureEJshows the angular dependence of the switching field 
for the full micromagnetic model for periods of 15, 25, 
50, and 100 ns and maximum applied field of 5 kOe at 
T = OK, and for periods of 15ns and 25ns with a max- 
imum applied field of lOkOe at T = 100 K. At OK, the 
general trend is for longer periods to reduce the switch- 
ing field. However, at 15°, the 100ns loop is observed to 
switch at a lower field than either the 50 or 25 ns loops. 
Similarly, at 30°, the 50 ns loop switches at a lower field 
than the 25 ns loop. One reason for this may be resonance 
in the switching fields for these angles and periods. At 
T = 100 K, the 25 ns loops switch at a lower field than the 
15 ns loops for all angles. At 90°, where thermal fluctu- 
ations are most prominent, the field at which relaxation 
occurs is independent of the period within the accuracy 
of the simulation. 

The increase of H sw with the misalignment angle in 
the micromagnetic simulation is consistent with recent 
experimental observations of Fe nanopillarsiSi 41 - 4 ^ 4 ^ 4 ^. 
However, the most recent experiment 4 ^ shows that a 
nanopillar with lateral dimension d ~ 5.2 nm, which 
our formulation suggests should show a dependence of 
_ff sw on misalignment angle similiar to the stack-of-spins 
model (i.e. like coherent rotation), actually exhibits the 
increasing dependence found in the full micromagnetic 
model. In addition, a nanopillar with lateral dimension 
d ~ 10 — 15 nm showed evidence of a multi-domain re- 
manence state. As noted in Ref. 44|> imperfections of the 



nanopillar structure appear to contribute to localized nu- 
cleation processes down to smaller than expected lateral 
dimensions, and probably also provide the pinning sites 
causing the multi-domain remanence state. This illus- 
trates the importance of coordinating experimental and 
simulation results in the micromagnetic approach. Fur- 
ther improvement of the predictions of the micromagnetic 
approach will likely have to incorporate such structural 
imperfections. 



IV. RECENT RESULTS FOR THE 2D KINETIC 
ISING MODEL 

Monte Carlo simulation of the Ising model, as well as 
other magnetic systems, continues to be an active field of 
research. Here, we present three recent results that are 
of interest in understanding the process of magnetization 
reversal in ultra-thin films. 



A. Dynamic Phase Transitions 

When the half-period ti/ 2 of the applied field is longer 
than the characteristic switching time in a constant field, 
(t(Hq)}, where Hq is the amplitude of the oscillating 
field, the magnetization can follow the changing field, re- 
sulting in standard hysteresis loops, such as those shown 
in Fig. n]and in Fig. However, when t\/% <C (r(iJ )), 
the magnetization cannot follow the field, but rather 
oscillates around one or the other of its zero-field sta- 
ble values. This breaking of the symmetry of the hys- 



(a) Snapshots of reversal at times (left to right) t = 7.350, 
7.375, 7.400, 7.425, and 7.450 ns 



(b) Snapshots of reversal at times (left to right) t = 0.00, 
3.000, 4.000, 6.250, and 8.500 ns 



FIG. 3: The z-component of the magnetization in the full micromagnetic model for the (a) 45° and (b) 90° hysteresis loops of 
Fig. la 



teresis loop is associated with a dynamic phase tran- 
sition (DPT) located at an intermediate value of the 
half-period. In terms of the dimensionless half-period, 
<9 = ti/2/ ( T (Ho)) , the transition is located at C rj 1. 
The dynamic order parameter for this transition is the 
period-averaged magnetization, 

2 r«( 2 *i/2) 
Qn = 7^7— / ro(t)dt . (13) 

-5*1/2 J(n-l)(2t 1/2 ) 

In Fig. we show hysteresis loops for the two- 
dimensional kinetic Ising model using Glauber dynamics 
for the dynamically disordered phase with 3> C and 
the dynamically ordered phase with 0C0 C . Time se- 
ries of Q n for 3> C , ~ 0c, and <C C are shown 
in Fig. E| 

The DPT was first discovered in numerical solutions 
of a mean-field model of a ferromagnet in an oscillating 
g e i c j45 j 46^ nag smce been intensively studied in mean- 
field modelo 4 7i 48 i 49 , kinetic Ising models 5n i 51 i 52 i 53 i 54 i 55 , 
the kinetic spherical model—, and anisotropic xySffi 
and HeisenbergS2s£R models. There have also been in- 
dications of its presence in experimental studies of hys- 
teresis in ultra-thin films of Cu on Co(001)£i*§£. From a 
theoretical point of view, its most interesting feature is 
that this far- from- equilibrium phase transition is a gen- 



uine continuous (second-order) phase transition that be- 
longs to the same universality class as the equilibrium 
phase transition in the Ising model in zero nc hi 53 - 54 - 55 - 63 . 
Unequivocal experimental verification of this interesting 
non-equilibrium phase transition is highly desirable and, 
given that new high-density magnetic recording media 
will require shorter reversal periods, may be relevant to 
the design of magnetic storage devices. 

B. Hysteresis Loop Area and Stochastic Resonance 

The values of Q measured for the stack-of-spins model 
described above appear to be consistent with the exis- 
tence of a DPT, although no detailed analysis has yet 
been made^. At lower frequencies, another interesting 
behavior is seen in both the stack-of-spins and kinetic 
Ising modelsS 4 ^. The normalized average hysteresis- 
loop area, 

is a measure of the average energy dissipation per period 
and is therefore a very important quantity. It is shown 
vs scaled frequency, 1/0, in Fig. OH for the stack-of-spins 
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FIG. 4: Angular dependence of the switching field for three models of magnetization reversal. The full micromagnetic model 
(diamonds) at T = 100 K shows a distinctly different behavior from both the stack-of-spins model at T — 20 K (circles) and 
the Stoner-Wohlfarth model (dashed line) 
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(a) Applied field periods are 15 ns (circles), 25ns (squares), 
50 ns (diamonds), and 100 ns (triangles). The field is sinusoidal 
with a maximum value of 5 kOe 



(b) Applied field periods are 15 ns (squares) and 25 ns (circles). 
The field is sinusoidal with a maximum value of 10 kOe 



FIG. 5: Angular dependence of the switching field for the full micromagnetic model at (a) OK and (b) 100 K. At OK, the LLG 
equation is completely deterministic, while at 100 K, it includes random fluctuations through a stochastic thermal field. 



model at T = 100 K and T = 20 K. At extremely low fre- 
quencies, the magnetization switches at very small values 
of H, so that (A) w 0. At high frequencies, the switch- 
ing rarely completes because the system is metastable for 
only a very short time interval. Therefore, M is nearly 
constant and again {A) w 0. A maximum in (^4) occurs 
at intermediate frequencies I/O w 0.1. For studies of 



hysteresis in a kinetic Ising model which switches by a 
single-droplet mechanism, this maximum was found to 
correspond to stochastic energy resonance^. This phe- 
nomenon has been studied further in the kinetic Ising 
mo ^ e r55 i 66 i 67^ anc j a i g0 recen tly investigated in models 

of superparamagnetic nanoparticles§& and Preisach sys- 
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FIG. 6: Simulated hysteresis loops for a kinetic Ising model (a) in the dynamically disordered phase for ^> C and (b) in 
the dynamically ordered phase for <C C . Data courtesy of S.W. Sides 
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First-Order Reversal Curves 



The First-Order Reversal Curve (FORC) technique 
was developed by Pike, et ali^ . in order to extract more 
information from magnetic samples than is represented 
by, for example, the coercive field or the remanent mag- 
netization. The FORC method has since been applied 
to a wide variety of systems, including several relevant 



to magnetic nanostructures£2iZ4iZ£i2iiZLZ§iZsL§P.. In addi- 
tion, progress has been made in understanding the role of 
reversible magnetization in the FORC method^! and in 
improving the efficiency of its computational use^ . Here, 
we illustrate the basic approach with an application to 
the kinetic Ising model. 

The FORC technique involves decreasing the applied 
field from a positive saturating field, Hq, to a series of 
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FIG. 8: Average hysteresis-loop area, (A), vs scaled frequency, 1/(9 for the stack-of-spins micromagnetic model. The same 
behavior is seen in two-dimensional Ising models that switch by a single-droplet mechanism, and the maximum is associated 
with stochastic energy resonance 



progressively more negative return fields, if r , and record- 
ing the normalized magnetization, m=M/Mg, as the 
field is increased from each of these return fields back 
to the positive saturating field. This process results in 
a family of first-order reversal curves, m(H r , H), where 
if represents the applied magnetic field as it is increased 
from H r back to Hq- Since the first-order reversal curves 
(FORCs) are determined by the type of reversal that has 
taken place before reaching if r , the full family of FORCs 
should contain useful information about the mechanisms 
of reversal. 

We can use the FORC method to better understand 
the process of hysteresis in the two-dimensional ferromag- 
netic kinetic Ising model on a square lattice, choosing the 
Glauber acceptance rule to produce the dynamic of the 
system with the energy given by (10). While most FORC 
studies have been done on systems with strong disorder, 
we focus here on the square-lattice Ising model without 
disorder. Our simulations were performed at a temper- 
ature of T = 0.8 T c which, given that k B T c « 2.269J 
for the two-dimensional square-lattice Ising model, cor- 
responds to ksT « 1.815 J. It has been foundS that the 
switching of a fully magnetized lattice for these param- 
eters occurs through single-droplet nucleation for fields 
up to | if | « 0.35, by multi-droplet nucleation for fields 
|if| Ri 0.35-0.9, and by strong-field (single-spin) rever- 
sal for fields |if | > 0.9. Since the process of switching 
is also influenced by the lattice size for finite lattices, 



these values serve only as guidelines. Here, we are mainly 
concerned with the multi-droplet regime, and so choose 
H = 0.55. 

We performed MC simulations to calculate the charac- 
teristic switching time r (for switching from m = 1.0 to 
to < —0.8) in a field of magnitude ifo = —0.55, finding 
t w 100 MCSS for a 128 x 128 lattice. We therefore 
chose a field period of P = 1000 MCSS, corresponding to 



a dimensionless half-period & 



P/2 



5. The form of 



the field is taken as a sawtooth, piecewise linear function 



H(t) = H 



i\t-P/2\ 
P 



(15) 



Figure shows the results of the simulation on a 128 
x 128 lattice for dimensionless half-periods of = 5, 10, 
and 25. The simulations were performed in parallel with 
100 independent realizations distributed over 20 proces- 
sors using the 48-bit linear congruential random number 
generator included with the SPRNG 2.0 package^. 

As the lattice just completes its reversal during the full 
hysteresis loop, we expect that the family of FORCs will 
reflect much of the dynamics that are occurring during 
the reversal. To investigate this, we divided the interval 
from if = [—0.55, 0.55] into 100 equal intervals. We be- 
gan the first FORC at a return field of if,, = 0.0, and 
recorded the magnetization at H a values corresponding 
to the endpoints of the 100 intervals. (Thus, for the first 
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FIG. 9: (a)Hysteresis loops for 0=5 (solid line), 10 (dashed line), and 25 (dotted line) on a 128 x 128 Ising lattice, (b) Family 
of FORCs for the same lattice with = 5 



FORC, wc took 51 values of the magnetization.) We then 
took a series of FORCs for H r values at the interval end- 
points from H = 0.0 to H = —0.55, producing a total 
of 51 FORCs. For each curve, we averaged over 100 re- 
alizations of the MC simulation, a technique commonly 
used to find the thermally averaged behavior of a system. 
The resulting family of FORCs is shown in Fig. 03. An 
animation of the reversal process for the FORCs shows 
that the reversal does proceed by the nucleation, growth, 
and shrinkage of multiple droplets (i.e., areas of reversed 
magnetization) . 

In a recent article^, we have continued this inves- 
tigation of the kinetic Ising model using the family of 
FORCs, as well as the FORC distribution, which can be 
derived from the FORCs as described in Ref. |7Q. The 
analysis yielded insights into the limits of application of 
the Kolmogorov-Johnson-Mehl-Avrami (KJMA) model 
of phase transformationSii2LS to the kinetic Ising sys- 
tem. In general, the FORC method appears to be quite 
sensitive to details of the magnetization reversal process, 
and with some thought can be helpful in developing in- 
sights into the construction of useful models. 



V. CONCLUSION 

Information storage devices utilizing magnetic nanos- 
tructures have become a technologically important part 
of our society. As demands for information storage in- 
crease, the size of the nanostructures must be decreased. 
At the same time, it becomes important to read and write 
the information to these devices (i.e. reverse the magne- 
tization) faster. The understanding of hysteresis in the 
magnetic nanostructures is therefore important to the 



continued growth of the information-storage industry. At 
the same time, the growth of computational resources has 
provided researchers with an invaluable tool with which 
to better understand these systems. 

In this overview, various common models and methods 
for simulating hysteresis in magnetic nanostructures have 
been presented along with results illustrating some of the 
properties of these systems. Micromagnetic simulations 
are accomplished by integration of the Landau-Lifshitz- 
Gilbert (LLG) equation. The LLG equation, despite be- 
ing both classical and phcnomcnological in origin, never- 
theless provides good insight into the magnetization dy- 
namics at nanosecond time scales, provided the system is 
sufficiently finely discretized. Our simulations on single 
Fe nanopillars show that the switching field (i.e. the field 
required to reduce Mz to 0) increases continuously as the 
angle between the z-axis and the applied field direction is 
increased, consistent with experiment. Reversal in these 
pillars is shown to nucleate at the endcaps and proceed 
by domain growth towards the center of the particle. The 
exception to this is the case of the applied field perpen- 
dicular to the long axis of the pillar, in which nucleation 
of reversal occurs along the whole length of the particle. 

Unfortunately, limitations on computer resources pre- 
vent extension of micromagnetic simulations beyond 
timescales of a few tens of nanoseconds. For timescales 
where the transition time for an individual spin to relax 
from the metastable to the stable state is much shorter 
than the time scale of interest, individual spin reversals 
occur with a probability which is related to the Boltz- 
mann factor. The dynamics of the system can then be 
modeled using kinetic Monte Carlo techniques with ei- 
ther the Ising or Heisenberg models. Here, we have shown 
three interesting applications of kinetic Monte Carlo sim- 
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ulations of a 2-D Ising model to understanding hystere- 
sis: dynamic phase transitions, stochastic resonance in 
the hysteresis loop area, and First-Order Reversal Curves 
(FORCs). These illustrate only a few of the ways sim- 
ulations of magnetic nanostructures may help give new 
insight into this important class of materials for ultra- 
high-density data storage. 
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